KSPZV1 correlate module expression, immunophenotypes, and cytokines to time to parasitemia and post-vax anti-CSP IgG

Objective:

To assess monotonic relationships between transcriptional modules, flow cytometric immunophenotypes, cytokines, CSP-specific IgG, and time-to-first parasitemia (up to 6 months).

Steps: 1. Filter and keep select features from non-transcriptomic data (FACS data and cytokine levels) 2. Load the expression set and filter the phenoData for Time point (0 for baseline, 25 for post-vaccination), CSP-specific IgG at 2 weeks post-vax, and time to parasitemia at 6 months. 3. Merge the selected phenoData features with the filtered non-transcritomic data. We will call this “alldata” 4. Collapse transcriptomic data (eset) into blood transcription modules (BTMs). 5. Merge the “alldata” with the BTMs or Monaco modules using patient ID. 6. For each dose group, build separate modules (one for each feature) adjusting for study site, presence of Pf infection at first vaccination, and pre-vax CSP-specific IgG 7. For each dose group, combine significant features in a unifying multi-variable model. 8. Convert pvalues to FDR and use it to keep only significant correlations

Notes:

Plasma cytokine concentrations were only included in the baseline but not in the correlations of post-vaccination and delta.

Load required packages

library(tidyverse)
library(tibble)
library(viridis)
library(ggpubr)
library(readxl)
library(ggcorrplot)
library(Biobase)
library(miscTools)
library(circlize)
library(RColorBrewer)
library(googledrive)
library(Hmisc)
library(tmod)

Import the expression set and non-transcriptomic data (FACS, cytokine, ICS)

##load full time to parasitemia data to 6 months for full 260 infants google drive
temp <- tempfile(fileext = ".rds")
dl <- drive_download(
  as_id("1OwUUrQZ_1jl_JExe3QWxCp-GaR83TtGG"), path = temp, overwrite = TRUE)
tte_data <- readRDS(file = dl$local_path)

##import all non-transcriptomic data google drive
temp <- tempfile(fileext = ".rds")
dl <- drive_download(
  as_id("1acv8D-gex3ps8o9SCWZT5ueP7XuPeHs3"), path = temp, overwrite = TRUE)
alldata <- readRDS(file = dl$local_path)

#load the baseline expression dataset
temp <- tempfile(fileext = ".rds")
dl <- drive_download(
  as_id("1togaBlNIxiDj16FyXTC-r7Qw0A9cG2az"), path = temp, overwrite = TRUE)
baseline_eset <- readRDS(file = dl$local_path)

#load the post-vax expression dataset
temp <- tempfile(fileext = ".rds")
dl <- drive_download(
  as_id("1w5Wzta8XISNHZ6KL-x9_qO98m3LxaOBa"), path = temp, overwrite = TRUE)
postvax_eset <- readRDS(file = dl$local_path)

#load the delta expression dataset
temp <- tempfile(fileext = ".rds")
dl <- drive_download(
  as_id("17xx0YaggLiyWdJc9rqA1nkPSE7dj05p0"), path = temp, overwrite = TRUE)
delta_eset <- readRDS(file = dl$local_path)

Convert gene expression data to BTM expression scores and join to other data

source("https://raw.githubusercontent.com/TranLab/ModuleLists/main/Gene2ModuleExpressionScores.R")

summary_stat <- median

#baseline
baseline_pas_hibtm <- Gene2ModuleExpressionScores(baseline_eset, module_list = "highBTMs", summary_stat = summary_stat)
rownames(baseline_pas_hibtm) <- paste(rownames(baseline_pas_hibtm), "highBTM_baseline")
rownames(baseline_pas_hibtm) <- gsub(" ", "_", rownames(baseline_pas_hibtm))

baseline_pas_monaco <- Gene2ModuleExpressionScores(baseline_eset, module_list = "MonacoModules", summary_stat = summary_stat)
rownames(baseline_pas_monaco) <- paste(rownames(baseline_pas_monaco), "monaco_baseline")
rownames(baseline_pas_monaco) <- gsub(" ", "_", rownames(baseline_pas_monaco))

baseline_pas_lowbtm <- Gene2ModuleExpressionScores(baseline_eset, module_list = "lowBTMs", summary_stat = summary_stat)
rownames(baseline_pas_lowbtm) <- paste(rownames(baseline_pas_lowbtm), "lowBTM_baseline")
rownames(baseline_pas_lowbtm) <- gsub(" ", "_", rownames(baseline_pas_lowbtm))

#post-vax
postvax_pas_hibtm <- Gene2ModuleExpressionScores(postvax_eset, module_list = "highBTMs", summary_stat = summary_stat)
rownames(postvax_pas_hibtm) <- paste(rownames(postvax_pas_hibtm), "highBTM_postvax")
rownames(postvax_pas_hibtm) <- gsub(" ", "_", rownames(postvax_pas_hibtm))

postvax_pas_monaco <- Gene2ModuleExpressionScores(postvax_eset, module_list = "MonacoModules", summary_stat = summary_stat)
rownames(postvax_pas_monaco) <- paste(rownames(postvax_pas_monaco), "monaco_postvax")
rownames(postvax_pas_monaco) <- gsub(" ", "_", rownames(postvax_pas_monaco))

postvax_pas_lowbtm <- Gene2ModuleExpressionScores(postvax_eset, module_list = "lowBTMs", summary_stat = summary_stat)
rownames(postvax_pas_lowbtm) <- paste(rownames(postvax_pas_lowbtm), "lowBTM_postvax")
rownames(postvax_pas_lowbtm) <- gsub(" ", "_", rownames(postvax_pas_lowbtm))

#delta
delta_pas_hibtm <- Gene2ModuleExpressionScores(delta_eset, module_list = "highBTMs", summary_stat = summary_stat)
rownames(delta_pas_hibtm) <- paste(rownames(delta_pas_hibtm), "highBTM_delta")
rownames(delta_pas_hibtm) <- gsub(" ", "_", rownames(delta_pas_hibtm))

delta_pas_monaco <- Gene2ModuleExpressionScores(delta_eset, module_list = "MonacoModules", summary_stat = summary_stat)
rownames(delta_pas_monaco) <- paste(rownames(delta_pas_monaco), "monaco_delta")
rownames(delta_pas_monaco) <- gsub(" ", "_", rownames(delta_pas_monaco))

delta_pas_lowbtm <- Gene2ModuleExpressionScores(delta_eset, module_list = "lowBTMs", summary_stat = summary_stat)
rownames(delta_pas_lowbtm) <- paste(rownames(delta_pas_lowbtm), "lowBTM_delta")
rownames(delta_pas_lowbtm) <- gsub(" ", "_", rownames(delta_pas_lowbtm))

#test innate only
data(tmod)
my_eset <- tmod # tmod includes Li et al. blood transcription modules and Chaussabel blood transcription modules.
baseline_eigengene_LI_DC <- eigengene(x = exprs(baseline_eset), g = fData(baseline_eset)$GeneSymbol, mset = my_eset, k = 1)
rownames(baseline_eigengene_LI_DC) <- paste(rownames(baseline_eigengene_LI_DC), "baseline", sep="_")
postvax_eigengene_LI_DC <- eigengene(x = exprs(postvax_eset), g = fData(postvax_eset)$GeneSymbol, mset = my_eset, k = 1)
rownames(postvax_eigengene_LI_DC) <- paste(rownames(postvax_eigengene_LI_DC), "postvax", sep="_")
delta_eigengene_LI_DC <- eigengene(x = exprs(delta_eset), g = fData(delta_eset)$GeneSymbol, mset = my_eset, k = 1)
rownames(delta_eigengene_LI_DC) <- paste(rownames(delta_eigengene_LI_DC), "delta", sep="_")

#remove _0 or _25 suffix from columnames
baseline_pas_all <- rbind(baseline_pas_hibtm, baseline_pas_monaco, baseline_pas_lowbtm, baseline_eigengene_LI_DC)
postvax_pas_all <- rbind(postvax_pas_hibtm, postvax_pas_monaco, postvax_pas_lowbtm, postvax_eigengene_LI_DC)
delta_pas_all <- rbind(delta_pas_hibtm, delta_pas_monaco, delta_pas_lowbtm, delta_eigengene_LI_DC)

colnames(baseline_pas_all) <- gsub("_.*", "", colnames(baseline_pas_all))
colnames(postvax_pas_all) <- gsub("_.*", "", colnames(postvax_pas_all))
colnames(delta_pas_all) <- gsub("_.*", "", colnames(delta_pas_all))

pas_all <- t(bind_rows(baseline_pas_all, postvax_pas_all, delta_pas_all)) %>%
  as.data.frame() %>%
  rownames_to_column(var = "PATID")

clin_dat <- pData(baseline_eset) %>%
  dplyr::select(PATID, treat, site, mal.vax.1, pfcsp_pre, pfcsp_post, log2FC_CSPAb)

alldata_select <- alldata %>%
  dplyr::select(PATID.OG, Timepoint, contains("FACS"), contains("CytokineObsConc")) %>%
  dplyr::rename(PATID = PATID.OG) %>%
  mutate(PATID = gsub("\\_.*", "", PATID)) %>%
  pivot_wider(., id_cols = PATID, names_from = Timepoint, values_from = 3:ncol(.),names_sep = "_") %>%
  select_if(~sum(!is.na(.)) > 0) %>% #remove columns with all NAs
  # dplyr::select(PATID, contains("FACS_PfSPZ"), matches("(FACS.*25)"), matches("(FACS_PfSPZ-specific_CD3.*25)"),
  #               matches("(FACS_CSP-specific.*25)"), contains("CytokineObsConc")) %>%
  dplyr::select(PATID, contains("FACS_"), contains("CytokineObsConc")) %>%
  dplyr::select(-c(contains("_of_TCRgd"))) %>% #remove redundant features (also represented by Vg9+Vd2+ of TCRgd)
  dplyr::select(-c(contains("of_live_lymphocytes"))) %>% #remove redundant features (also represented by live_leukocytes)
  dplyr::select(-c(contains("of_live_monocytes"))) %>% #remove redundant features (also represented by live_leukocytes)
  dplyr::select(-c(contains("CSP-spec_"))) #remove redundant features (also represented by live_leukocytes)

all_data_btm <- tte_data %>%
  full_join(., clin_dat, by = "PATID") %>%
  full_join(., alldata_select,  by = "PATID") %>%
  full_join(., pas_all,  by = "PATID") %>%
  mutate(treat = factor(treat, levels = c("Placebo", "4.5 x 10^5 PfSPZ", "9.0 x 10^5 PfSPZ", "1.8 x 10^6 PfSPZ"))) %>%
  mutate(mal.vax.1 = factor(ifelse(mal.vax.1 == 0, "pfneg", "pfpos"))) %>%
  mutate(site = factor(site))  %>%
  dplyr::select(PATID, treat, site, tte.mal.atp.6, everything())

The following is applicable to baseline data only: columns that have a large percentage of zero values and a few high leverage points tend to show up as significant correlations downstream. To avoid this, remove columns that contain ~>66% zero values. Only certain cytokines fell into this category, and by thresholding, we were able to remove them and prevent ‘false’ correlations.

# not needed for spearman correlation
ncol(all_data_btm)

## [1] 3046

i <- colSums(all_data_btm == 0, na.rm=TRUE) < 0.66*nrow(all_data_btm)
all_data_btm <- all_data_btm[i]
ncol(all_data_btm)

## [1] 3040

Arrange module expression data for pairwise correlations and linear models

Low-annotation BTMs collapsed by median of all members chose in manuscript for all analyses. This method has been used in other systems vaccinology studies over other approaches such as eigenvalues (first principal component) or other measures of variance.

Note here that we include baseline plasma cytokine data with the module expression as this was not included in any prior analyses.

#reduce to low annotation BTMs
#LI. and DC. modules are eigenvalue-collapsed modules from tmod
#_monoaco, _highBTMs, and _lowBTMs are median-collapsed modules

all_data_btm_reduced <- all_data_btm %>%
  dplyr::select(-c(contains("DC."), contains("_highBTM"), contains("LI."), contains("_monaco"))) #remove other modules

all_module_lm_dat <- all_data_btm_reduced %>%
  dplyr::select(PATID, treat, site, mal.vax.1, tte.mal.atp.6, everything()) %>%
  pivot_longer(., cols = 6:ncol(.), names_to = "feature", values_to = "value") %>%
  mutate(feature = gsub("_0", "_baseline", feature)) %>% #recode 0 as baseline
  mutate(feature = gsub("_25", "_postvax", feature)) %>% #recode 25 as postvax
  mutate(feature = gsub("pfcsp_pre", "CSP-specific IgG_baseline", feature)) %>%
  mutate(feature = gsub("pfcsp_post", "CSP-specific IgG_postvax", feature)) %>%
  mutate(feature = gsub("log2FC_CSPAb", "CSP-specific IgG_delta", feature)) %>%
  mutate(timepoint = "TBD") %>% #features are TBD until specified
  mutate(timepoint = ifelse(grepl("baseline", feature), "baseline", timepoint)) %>% 
  mutate(timepoint = ifelse(grepl("postvax", feature), "postvax", timepoint)) %>%
  mutate(timepoint = ifelse(grepl("delta", feature), "delta", timepoint)) %>%
  mutate(feature = gsub("_baseline", "", feature)) %>% #remove suffix
  mutate(feature = gsub("_postvax", "", feature)) %>%
  mutate(feature = gsub("_delta", "", feature)) %>%
  mutate(feature = gsub("_lowBTM", "", feature)) %>%
  mutate(feature = gsub("_highBTM", "", feature)) %>%
  mutate(feature = gsub("_monaco", "", feature)) %>%
  pivot_wider(., names_from = timepoint, values_from = value) %>% #pivot to wide format
  dplyr::select(PATID, treat, site, mal.vax.1, tte.mal.atp.6, feature, baseline, postvax, delta) %>%
  mutate(delta = ifelse(grepl("FACS", feature), log2((postvax+1e-06)/(baseline+1e-06)), delta))

Perform pairwise correlations between each feature and time to parasitemia and fold-change CSP-specific IgG, all dose groups combined

library(Hmisc)
my_significance_thres <- 0.05
flattenCorrMatrix <- function(nmat, cormat, pmat) {
  ut <- upper.tri(cormat)
  data.frame(
    row = rownames(cormat)[row(cormat)[ut]],
    column = rownames(cormat)[col(cormat)[ut]],
    n = nmat[ut],
    cor  =(cormat)[ut],
    p = pmat[ut]
    )
}


eigen_btm_cor_dat_all <- all_module_lm_dat %>%
  dplyr::select(PATID, treat, tte.mal.atp.6, feature, baseline, postvax, delta) %>%
  pivot_longer(., cols = c(baseline, postvax, delta), names_to = "timepoint", values_to = "value") %>%
  dplyr::select(PATID, timepoint, treat, tte.mal.atp.6, everything()) %>%
  mutate(feature = plyr::mapvalues(feature, from = tmod$gs$ID, to = paste0(tmod$gs$Title, " (",tmod$gs$original.ID, ")"), warn_missing =F))  %>%
  filter(!grepl("TBA", feature)) %>%
  drop_na(value, tte.mal.atp.6) %>%
  pivot_wider(., names_from = feature, values_from = value) %>%
  mutate(timepoint = factor(timepoint, levels = c("baseline", "postvax", "delta"))) %>%
  arrange(timepoint, treat, PATID)

eigen_btm_cor_dat_all_CSP_IgG <- eigen_btm_cor_dat_all %>%
  dplyr::select(PATID, timepoint, treat, "CSP-specific IgG") %>%
  pivot_wider(names_from = "timepoint", values_from = "CSP-specific IgG", names_prefix = "CSP_specific_IgG_")

eigen_btm_cor_dat_all <- eigen_btm_cor_dat_all %>%
  full_join(., eigen_btm_cor_dat_all_CSP_IgG, by = c("PATID", "treat")) %>%
  dplyr::select(PATID, timepoint, treat, tte.mal.atp.6, contains("CSP_specific_IgG"), everything()) %>%
  dplyr::select(-c("CSP-specific IgG", "CSP_specific_IgG_baseline", "CSP_specific_IgG_delta"))

temp_df <- corres_list <- flatten_cormat_list <- sig_features_list <- significant_cormat_r <- significant_cormat_p <- significant_cormat_fdr <- c()
significant_cormat_n <- significant_cormat_r_networkplot <- c()
for(i in unique(eigen_btm_cor_dat_all$timepoint)){
    temp_df <- eigen_btm_cor_dat_all %>%
      filter(timepoint == i) %>%
      dplyr::select(-c(timepoint, treat)) %>%
      column_to_rownames(var = "PATID") 
    corres_list[[i]]$n <- rcorr(as.matrix(temp_df), type = "spearman")$n
    corres_list[[i]]$r <- rcorr(as.matrix(temp_df), type = "spearman")$r
    corres_list[[i]]$p <- rcorr(as.matrix(temp_df), type = "spearman")$P
    corres_list[[i]]$fdr <- tabletools::rcorr_padjust(rcorr(as.matrix(temp_df), type = "spearman"), method = "BH")$P
    corres_list[[i]]$flattened_cormat <- flattenCorrMatrix(corres_list[[i]]$n,
                                                                corres_list[[i]]$r,
                                                                corres_list[[i]]$p)
    corres_list[[i]]$flattened_cormat$fdr <- p.adjust(corres_list[[i]]$flattened_cormat$p, method = "BH") #use BH for FDR correction
    corres_list[[i]]$flattened_cormat <- corres_list[[i]]$flattened_cormat[complete.cases(corres_list[[i]]$flattened_cormat),]
    flatten_cormat_list[[i]] <- corres_list[[i]]$flattened_cormat
    sig_features_list[[i]] <- corres_list[[i]]$flattened_cormat %>%
      filter(row == "tte.mal.atp.6" & fdr < my_significance_thres |
               column == "tte.mal.atp.6" & fdr < my_significance_thres |
               row == "CSP_specific_IgG_postvax" & fdr < my_significance_thres |
               column == "CSP_specific_IgG_postvax" & fdr < my_significance_thres) %>%
      dplyr::select(row, column, n, cor, p, fdr) %>%
      dplyr::rename(feature1 = "row") %>%
      dplyr::rename(feature2 = "column")
      significant_features <- unique(c(sig_features_list[[i]]$feature1, sig_features_list[[i]]$feature2))
      significant_cormat_r[[i]] <- corres_list[[i]]$r[significant_features, significant_features] #make downselected correlation matrix
      significant_cormat_n[[i]] <- corres_list[[i]]$n[significant_features, significant_features] #make downselected sample size matrix
      significant_cormat_p[[i]] <- corres_list[[i]]$p[significant_features, significant_features] #make downselected p value matrix
      significant_cormat_fdr[[i]] <-  corres_list[[i]]$fdr[significant_features, significant_features] #make downselected fdr matrix
      significant_cormat_r_networkplot[[i]] <- significant_cormat_r[[i]]
      significant_cormat_r_networkplot[[i]][significant_cormat_fdr[[i]] >= my_significance_thres] <- 0 #for FDR >= threshold, set to 0
    }

## Warning in sqrt(npair - 2): NaNs produced

## Warning in sqrt(npair - 2): NaNs produced

## Warning in sqrt(npair - 2): NaNs produced

## Warning in sqrt(npair - 2): NaNs produced

## Warning in sqrt(npair - 2): NaNs produced

## Warning in sqrt(npair - 2): NaNs produced

## Warning in sqrt(npair - 2): NaNs produced

## Warning in sqrt(npair - 2): NaNs produced

#bind correlation dataframes
all_cor_dfs <- bind_rows(flatten_cormat_list, .id = "timepoint") %>%
  dplyr::select(timepoint, everything()) %>%
  mutate(treat = factor(timepoint, levels = c("baseline", "postvax", "delta")))

#bind signifificant features
all_sig_features <- bind_rows(sig_features_list, .id = "timepoint") %>%
  dplyr::select(timepoint, everything()) %>%
  mutate(treat = factor(timepoint, levels = c("baseline", "postvax", "delta")))

Clean row and columns names to plot data as network plots

for(i in names(significant_cormat_r)){
  #r
  colnames(significant_cormat_r[[i]]) <- rownames(significant_cormat_r[[i]]) <- gsub("tte.mal.atp.6", "time to parasitemia",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r[[i]]) <- rownames(significant_cormat_r[[i]]) <- gsub("CSP_specific_IgG_", "CSP-specific IgG_",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r[[i]]) <- rownames(significant_cormat_r[[i]]) <- gsub("FACS_|_innate_BTM|_highBTM|_monaco|CytokineObsConc_", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r[[i]]) <- rownames(significant_cormat_r[[i]]) <- gsub("_of_live_PBMCs", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r[[i]]) <- rownames(significant_cormat_r[[i]]) <- gsub("_of_live_leukocytes", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r[[i]]) <- rownames(significant_cormat_r[[i]]) <- gsub("PfSPZ", "Pf",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r[[i]]) <- rownames(significant_cormat_r[[i]]) <- gsub("_", " ",
                                                                                     rownames(significant_cormat_r[[i]]))
  #network r
  colnames(significant_cormat_r_networkplot[[i]]) <- rownames(significant_cormat_r_networkplot[[i]]) <- gsub("tte.mal.atp.6", "time to parasitemia",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r_networkplot[[i]]) <- rownames(significant_cormat_r_networkplot[[i]]) <- gsub("CSP_specific_IgG_", "CSP-specific IgG_",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r_networkplot[[i]]) <- rownames(significant_cormat_r_networkplot[[i]]) <-
    gsub("FACS_|_innate_BTM|_highBTM|_monaco|CytokineObsConc_", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r_networkplot[[i]]) <- rownames(significant_cormat_r_networkplot[[i]]) <- gsub("_of_live_PBMCs", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r_networkplot[[i]]) <- rownames(significant_cormat_r_networkplot[[i]]) <- gsub("_of_live_leukocytes", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r_networkplot[[i]]) <- rownames(significant_cormat_r_networkplot[[i]]) <- gsub("PfSPZ", "Pf",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_r_networkplot[[i]]) <- rownames(significant_cormat_r_networkplot[[i]]) <- gsub("_", " ",
                                                                                     rownames(significant_cormat_r[[i]]))
  #fdr
  colnames(significant_cormat_fdr[[i]]) <- rownames(significant_cormat_fdr[[i]]) <- gsub("tte.mal.atp.6", "time to parasitemia",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_fdr[[i]]) <- rownames(significant_cormat_fdr[[i]]) <- gsub("CSP_specific_IgG_", "CSP-specific IgG_",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_fdr[[i]]) <- rownames(significant_cormat_fdr[[i]]) <- gsub("FACS_|_innate_BTM|_highBTM|_monaco|CytokineObsConc_", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_fdr[[i]]) <- rownames(significant_cormat_fdr[[i]]) <- gsub("_of_live_PBMCs", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_fdr[[i]]) <- rownames(significant_cormat_fdr[[i]]) <- gsub("_of_live_leukocytes", "",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_fdr[[i]]) <- rownames(significant_cormat_fdr[[i]]) <- gsub("PfSPZ", "Pf",
                                                                                     rownames(significant_cormat_r[[i]]))
  colnames(significant_cormat_fdr[[i]]) <- rownames(significant_cormat_fdr[[i]]) <- gsub("_", " ",
                                                                                     rownames(significant_cormat_r[[i]]))
  }

Plot correlation matrix, all dose groups

Only show features that significantly correlate with either time to parasitemia or log-fold-change CSP-specific IgG

library(corrr)
min_cor <- 0.1
tl.cex <- 9
lab_size <- 1.65

my_signif_cormat_alldose_baseline <- ggcorrplot(corr = significant_cormat_r$baseline,
                                                method = "square",
                                                hc.order = FALSE,
                                                hc.method = "ward.D2",
                                                outline.col = "white",
                                                ggtheme = ggplot2::theme_minimal(),
                                                colors = c("darkblue", "white", "darkred"),
                                                type = "lower",
                                                sig.level = my_significance_thres,
                                                p.mat = significant_cormat_fdr$baseline,
                                                insig = "blank",
                                                lab = "true",
                                                tl.cex = tl.cex,
                                                lab_size = lab_size)

my_signif_networkplot_alldose_baseline <- network_plot(significant_cormat_r_networkplot$baseline, 
                                legend = "full",
                                colors = c("darkblue", "white", "darkred"),
                                min_cor = .1,
                                curved = TRUE)

my_signif_cormat_alldose_postvax <- ggcorrplot(corr = significant_cormat_r$postvax,
                                                method = "square",
                                                hc.order = FALSE,
                                                hc.method = "ward.D2",
                                                outline.col = "white",
                                                ggtheme = ggplot2::theme_minimal(),
                                                colors = c("darkblue", "white", "darkred"),
                                                type = "lower",
                                                sig.level = my_significance_thres,
                                                p.mat = significant_cormat_fdr$postvax,
                                                insig = "blank",
                                                lab = "true",
                                                tl.cex = tl.cex,
                                                lab_size = lab_size)

my_signif_networkplot_alldose_postvax <- network_plot(significant_cormat_r_networkplot$postvax, 
                                legend = "full",
                                colors = c("darkblue", "white", "darkred"),
                                min_cor = .1,
                                curved = TRUE)


my_signif_cormat_alldose_delta <- ggcorrplot(corr = significant_cormat_r$delta,
                                                method = "square",
                                                hc.order = FALSE,
                                                hc.method = "ward.D2",
                                                outline.col = "white",
                                                ggtheme = ggplot2::theme_minimal(),
                                                colors = c("darkblue", "white", "darkred"),
                                                type = "lower",
                                                sig.level = my_significance_thres,
                                                p.mat = significant_cormat_fdr$delta,
                                                insig = "blank",
                                                lab = "true",
                                                tl.cex = tl.cex,
                                                lab_size = lab_size)

my_signif_networkplot_alldose_delta <- network_plot(significant_cormat_r_networkplot$delta, 
                                legend = "full",
                                colors = c("darkblue", "white", "darkred"),
                                min_cor = .1,
                                curved = TRUE)

Correlation matrices (corresponds to Figure S6A)

Spearman’s correlation, FDR<5%

ggarrange(my_signif_cormat_alldose_baseline, my_signif_cormat_alldose_postvax, my_signif_cormat_alldose_delta,
          common.legend = TRUE, align = "hv",
          labels = c("baseline","2 weeks post-vax", "delta"),
          nrow=3,
          heights = c(3,4,4))

BTM-level network plots (Figure 5A of manuscript)

Spearman’s correlation, FDR<5%

ggarrange(my_signif_networkplot_alldose_baseline, my_signif_networkplot_alldose_postvax, my_signif_networkplot_alldose_delta,
          common.legend = TRUE, align = "hv",
          labels = c("baseline","2 weeks post-vax", "delta"),
          nrow=3,
          widths = c(1,1.27,1.27))